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We present a new numerical algorithm for the solution of coupled col- 
lisional and collisionless systems, based on the block structured adaptive 
mesh and time refinement strategy (AMR). We describe the issues associ- 
ated with the discretization of the system equations and the synchroniza- 
tion of the numerical solution on the hierarchy of grid levels. We implement 
a code based on a higher order, conservative and directionally unsplit Go- 
dunov's method for hydrodynamics; a symmetric, time centered modified 
symplectic scheme for collisionless component; and a multilevel, multigrid 
relaxation algorithm for the elliptic equation coupling the two components. 
Numerical results that illustrate the accuracy of the code and the relative 
merit of various implemented schemes are also presented. 
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1. INTRODUCTION 

Astro-physical systems are typically complex and highly nonlinear, providing the 
ground for the occurrence, either singly or concomitantly, of numerous physical pro- 
cesses. These include, among others, hydro- and magnetohydro-dynamics, gravity, 
radiation and many-body interactions. They operate on a wide range of spatial and 
temporal scales and it is often desirable to fully cover these ranges for a thorough 
understanding of the problem. Thus the problems are demanding both in terms 
of physics algorithms and dynamic range (resolution). While the development of 
high order numerical schemes certainly improves the quality of numerical solutions 
and the availability of ever more powerful computers has allowed the performance 
of larger calculations, special techniques are required in order to achieve very large 
dynamic ranges. 
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The Adaptive Mesh Refinement (AMR) technique offers a powerful solution for 
this purpose [1, 2]. While there are difficulties associated with its implementation, 
the use of AMR in astrophysics and cosmology has grown significantly to include 
studies of nucleosynthesis in Supernovae explosions [3, 4, 5, 6], of the multiphase 
interstellar medium and radiative shock hydrodynamics [7, 8] , the problem of star 
formation out of the collapse of protostellar clouds [9, 10, 11] and the formation 
of the first stars as well as the large scale structure in the universe [12, 13]. The 
use of AMR technique in the above examples has often been instrumental in either 
revealing new properties of the investigated system (such as instabilities) or pointing 
out mistaken views based on limited resolution calculations. 

We will consider two aspects of application of AMR to astrophysical problems. 
The first is the extension to incorporate self-gravity. This capability leads to new 
algorithmic difficulties due to the elliptic nature of the problem. In particular 
the solution has to be computed simultaneously on all levels of refinement and 
continuity of both the solution and its normal derivative has to be enforced at 
coarse- fine level interfaces [14]. However, this introduces non-trivial complications 
when time refinement is also employed: the coarser levels are advanced first and 
finer levels are advanced with the assumption that boundary conditions at fine- 
coarse interfaces are provided by the coarse level solutions interpolated in time. 
However, the full multilevel elliptic solution can only be computed when all levels 
are synchronized and thus is not available when the coarser levels are ahead of the 
finer levels. Thus the first implementations of a full multilevel elliptic solver for 
self-gravity [10, 11, 15] do not use refinement in time. And when employing time 
refinement the multilevel elliptic equation has been solved as a set of independent 
boundary value problems, one for each level, not a fully multilevel solution [16, 17, 
18, 19]. 

The second issue we will address is the application of AMR to hybrid systems, 
that is a self-gravitating gas coupled to a particle representation of a collisionless 
matter described by Vlasov-Poisson equations [20] . This is relevant to several 
problems in astrophysics, particularly for modeling the formation and evolution 
of structure in the universe. In this case the AMR technique is combined with 
Particle-Mesh methods to compute the right-hand side to the Poisson's equation 
due to the particles mass. In order to take advantage of the higher resolution of the 
finer grids, it is desirable to advance the particles with the force compute on, and 
the timesteps of, the finest grids that cover their spatial position. This introduces 
a complication for the elliptic solver similar to the one described above for the case 
of self-gravitating gas dynamics, in the sense that the particles contribution to the 
right-hand side of Poisson's equation needs to be accounted for even when the levels 
of refinement are not synchronized. 

Various AMR codes for hybrid systems have been developed over the past several 
years [21, 16, 17, 18, 19, 15]. The refinement strategy in [16, 18, 15] is based on 
splitting of individual cells and in [17] only the collisionless component is evolved. 
Virtually all schemes use Strang splitting [22] for the multidimensional version 
of the hydrodynamics and a modified leap-frog method for the integration of the 
equation of motion of the particles. This has also been done for problems in which 
there is no coupling to a collisional fluid component, as occurs in computations of 
collisionless plasmas [23]. 
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In our approach we use the block structured scheme for adaptive mesh and time 
refinement proposed in [2] as a starting point and extend it include gravity and 
collisionless particle dynamics. We use an unsplit Godunov's method for hydrody- 
namics [24]; a symmetric, time centered modified symplectic scheme based on the 
kick-drift-kick sequence for the collisionless component; and a multilevel, multigrid 
relaxation algorithm for the elliptic equation coupling the two systems. We in- 
troduce two new procedures to solve synchronization issues described above that 
arise with the elliptic solver when the coarse and finer levels are not synchronized. 
We use a method analogous to those developed for AMR for incompressible flows 
[25, 26] to compute a lagged estimate of the correction of the elliptic matching 
conditions at boundaries between refinement levels at times when the levels are not 
synchronized. We also present a detailed discussion of refinement in time in the 
presence of collisionless particles, including methods for associating particles with 
refinement levels, and a particle aggregation operation to cost-effectively estimate 
the density distribution on coarse levels due to particles evolved on finer levels 
without compromising the code accuracy. We provide a detailed description of the 
formal discretization of the system of equations and the issues involved with the 
synchronization of the numerical solution in the presence of refinement in time, at 
a level of detail which we feel is lacking in the current literature. 

The paper is organized as follows. First, the evolution equations and single level 
algorithms are outlined in Section 2. In Section 3 we provide a formal definition of 
the employed AMR volume discretization, variables and operators and then describe 
in detail the general AMR algorithm for hybrid system. Finally, in section 4 we 
test the accuracy of the code by comparing its results against for set of standard 
solutions. 

2. EVOLUTION EQUATIONS AND TEMPORAL 
DISCRETIZATION 

In this section we introduce the system equations and describe their temporal 
discretization on individual levels of refinement. Motivated by cosmological applica- 
tions, the numerical schemes are formulated for a grid with a time dependent scale 
length, a(t). Thus after a description of the expanding grid, we introduce the time 
discretization for the equations of hydrodynamics with gravity and the equations 
of motion for the collisionless component and then briefly outline the time step 
constraints and code units. 

2.1. Comoving Frame 

Cosmic expansion is described by the first Friedmann's equation which reads [27] 

- = H (O m a' 3 + n k a' 2 + A ) 1/2 (1) 

where a(t) is the scale of the universe as a function of time; H measures the 
current (a = 1) rate of cosmic expansion (the Hubble constant); fi m , fifc and Q\ are 
parameters representing the current energy density associated to matter, curvature 
and 'dark' component, respectively, in units of the closure value. The solution to 
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Eq. (1) relating cosmic time and expansion parameter reads 

#o t(a) = P m (2) 

Jo a (fi m a~ 3 + flk a~ 2 + S1a) 

and admits simple solutions for Qk = 0. Since the equations of motion in an 
expanding Universe are most naturally solved in a comoving frame, which expands 
at the rate a/ a given by Eq. (1), we operate the change of coordinates 



x 



aity 1 r (3) 



where r and x are the coordinates in the laboratory and comoving reference frame 
respectively, and transform all differential operators (time derivatives, gradient, 
laplacian) according to [27] 
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(4) 

d_\ f d_\ dx /d_\ _(d\ _axfd_ 



(5) 



The velocity 



a . , . 

r = ax + ax = —r + ax (b) 

is decomposed into a Hubble flow, (a/ a) r, and a peculiar proper component, u — 
ax. It is also convenient to introduce the density and pressure in terms of the 
comoving volume X j clS opposed to the proper volume r , 

p(t,x = r/a) = a 3 p p (t,r) (7) 
P(t,x = r/a) = a 3 P p (t,r) (8) 

where the subscript p indicates the proper quantities. 

2.2. Hydrodynamics 

In comoving coordinates, the hydrodynamics is described by the following set of 
inhomogeneous partial differential equations 

dpui Id. a 1 d(j) 

^T + ad^ k {pU * Uk + P5lk) - -a PUi -a P d^ (10) 

doe Id.. „. . „ct 1 dS 

+ Pe + PK = -2-pe--pu k -^ 11 

ot a dxk a a dxk 

dps Id,. „d 
dt a dxk a 

expressing (from top to bottom) mass, momentum, energy and entropy conser- 
vation. Here, p and P are the comoving density and pressure respectively, u is 
the peculiar proper velocity, <j) the proper gravitational potential, e = e t h + &k — 
P/p(l — 1) + m 2 /2 is the specific total energy, s = P/p 7 is the specific entropy and 



! Please write \titlerunninghead{< (Shortened) Article Title>} in file ! 



•5 



7 the gas adiabatic index. The first inhomogeneous terms on the right hand side 
of Eq. (10)-(12) describe the effects due to adiabatic cosmic expansion (cx a/ a). In 
particular, the factor 2 for the last two equations arises by assuming that the in- 
ternal energy of the gas is solely associated with translational degrees of freedom 1 . 
Finally, gravity is described by the gravitational potential <fi, which is generated by 
the matter distribution of both the collisional and the collisionless components. 

We use a cell-centered discretization for our primary dependent variables: U(x, t) = 
(p, pu, pe) T — ► f™, where i £ Z D indexes grid points in a space with D dimensions 
and n is a discrete time index. The starting point for our temporal discretization 
is a conservative finite-difference method for the hydrodynamic equations: 



U 



n+l,h 



= U n - AtA{t) (d ■ F n+ ^ (13) 



where D ■ F n+ 2 approximates the spatial derivative terms on the left-hand side of 

! ! „n + l/2 „n + l/2 



(11), A{t) is a diagonal matrix with elements a „+i/ 2 > »+i ; ( a ™+iyi 1 ^rrrp , an d 
At is the time step. We use an unsplit second-order Godunov's method [28, 29, 24] 
to compute the flux divergence. In addition to the cosmological expansion terms, 
the gas and particle components couple through the force field which is solution to 
the following Poisson's equation 

A<\> = —— (p m - ( Pm )) . (14) 

Here p m = p gas + p palt is the total comoving mass density; the particle density, 
Ppart, is computed through a Particle Mesh method (see details in the Appendix). 
When periodic boundary conditions are used, (p m ) is the volume average, otherwise 
it is zero. The details of the spatial discretization for Poisson's equation are given 
in the next section. For now, we assume that we can compute the gravitational 
acceleration at cell centers, f t w — V0, with second order accuracy. 

To compute the effect of the source terms, we need to compute them before and 
after the hydrodynamic update, taking advantage of the fact that no sources appear 
in the density evolution equation: 



S? = 



/ \ 

p n+l/2 f n /an+1 +p n u n [{a n /a n+l ) _ 1] 

A(pe fe ) + PP [{a n /a n+1 f - 1] 
V P?s? [(a"/« n+1 ) 2 - 1] J 



(15) 



where ^ +1/2 = (1/2) (p? + /^ +1 ), and A(pe k ) = (1/2) K +1 « +1 ) 2 - p^) 2 ]. 
After the hydrodynamics update given in Eq.(13) we apply the source terms as 
follows 

U n+1 ' h ' Sl = U n+1 ' h + AtS n . (16) 

The above estimate of the source term, after being converted to primitive variable 
form, is also used in the hydrodynamic predictor step in order to obtain fully 



1 More generally, in D dimensions, the energy losses due to expansion are eth/ e th = f (7—1)0/0 
for the internal energy and ej./ej. = 2d /a for the specific kinetic energy respectively. 
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time-centered fluxes. S n accounts for the expansion terms to the desired accuracy 
but is only first order accurate as far as the gravity term is concerned (hence the 
superscript Si). After the new gravitational potential has been computed a source 
correction term is estimated as 



SS, 



™+3 _ 1 „n+l 5 f 



n+1/ 2 



( 



2« 



,n+l 



+ [At/4 



\ 



i n+1 ] 5f" +1/2 



(17) 



with, 5fi +1 ^ = fl +1 — and the final source update is 



jjn+l = jjn+l.Mi + AtSS n+ ^. 



(18) 



2.2.1. Hypersonic Flows 

Accretion flows induced by gravity are typically hypersonic and can be char- 
acterized by very large Mach numbers M > 100. This situation is common in 
cosmological simulations [30] . 

In this case the total energy is largely dominated by the kinetic component, 
efe ~ M 2 eth- Since conservative hydro-schemes track the total energy, e = e t u + ek, 
relatively small errors in the partition of the two components can produce spurious 
values of e t h- This is particularly worrisome when a 4-byte (single precision) digital 
representation of the numerical data is employed. For this reason, we introduce 
the additional equation (12) describing the evolution of the gas entropy. When the 
Mach number of the bulk flow is very high, M > 50, and away from shocks, it 
provides a more accurate solution from the thermal energy than the total energy 
equation. Eq. (12) is naturally incorporated in the numerical scheme for hydrody- 
namics adopted here because of its conservative form. In addition, being a simple 
advection equation, the conservative fluxes for its integration are a byproduct of 
the Riemann solution and require virtually no extra effort to compute. 

The equation for the internal energy could be alternatively integrated [31], but 
its non-conservative form makes it less attractive. Authors in Ref. [32] already 
employed the entropy equation in order to improve the accuracy of their Total 
Variation Diminishing scheme for hypersonic flows, although their implementation 
required solving twice for the hydrodynamic equations with an extra cost of 30%. 

It is worth pointing out that for high Mach number flows, errors in the velocity 
field are propagated into the Riemann solver solution for the pressure at the cell 
interface, P*, with a coefficient ~ A4, that is, largely amplified. More precisely we 
find 

— = 7 — M (19) 

where Au is the one dimensional velocity jump at the cell interface. However, 
as illustrated above, the pressure terms enter the hydrodynamic equations with a 
weight M^ 1 as compared to kinetic terms and, therefore, the numerical solution 
is not degraded in the case of high Mach number flows because of approximations 
involved in the Riemann solver solution. 



! Please write \titlerunninghead{< (Shortened) Article Title>} in file ! 



7 



2.3. Collisionless Component 

The collisionless component is described by a set of particles whose evolution in 
phase space is computed according to 

d J= -u (20) 
dt a v ' 

du a 1 „ 

-77 = —U+-f 21 

at a a 

where y and u are the comoving coordinate and peculiar proper velocity respec- 
tively. The acceleration acting on the particle, /, is obtained by first computing 
the acceleration on the grid, using a cell- or face- centered scheme, and then by 
interpolating it to the particle position through a Particle Mesh method [20] . 

In order to advance in time the particle positions and velocities we propose the 
following integration scheme based on a kick-drift-kick sequence [33] : first the par- 
ticles velocity and positions are updated as 

«" +V2 = »"^ + ^ny")y (22) 



n+1 = y n + -^u n+1 / 2 At. (23) 



V y ' a n+l/2 

After computing the acceleration at the new timestep, the particle velocity is finally 
updated as 

= «" +1/2 ^? + Jti r +1 (v n+1 )^- (24) 



The proposed scheme is reflexive and hence symplectic [34]. This has the nice 
property that the integral of motion will be conserved on average preventing sec- 
ular accumulation of error and keeping the system about its true trajectory in 
phase space (see, e.g., discussion in [33]). We have also implemented an alterna- 
tive method, based on the more common drift-kick-drift sequence, which does the 
following: first particle positions at half the time step are predicted as 



y n+l/2 n + J_ »^ 

y y t- a n 2 ' 



(25) 



then particle positions and velocity are further temporarily updated to 

1 



u n+l.* = u n 



r n+l 



y n+l.* _ yn+1/2 



1 



■.n+1 



U 



n+1,* 



At 



(26) 
(27) 



Based on the gravitational potential at the new time-step, a final correction term 
is applied, that allows second order accuracy 



u 



n+1 



,n+l 



f n+1 (y n+1/2 )-f n (y n+1/2 ] 



At 
~2~ 



y n+l = y n+l/2 



n n+l 



U 



n+1 



At 



(28) 
(29) 
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The above scheme, however, is not fully reflexive. The need of temporary states 
that approximate the solution at the end of the time step in order to calculate the 
final correction step, breaks the time symmetry of the scheme. 

Note that in both schemes, there is no need for storage of extra information about 
either the particle positions or their velocities at any old or intermediate step. In 
addition, the gravitational potential is computed only once per time step making 
their overall computational cost of rather inexpensive. 

It should be noticed, however, that in the case of AMR reflexivity is lost even in 
the former scheme without additional precautions. This occurs when a particle is 
transferred to a different level of refinement, or even during a refinement operation. 
In both cases, the gravitational potential changes as a result of these procedures. 
Since this change takes place after the correction steps, the backward application 
of the scheme would not reproduce the initial configuration. Although in principle 
this could be fixed by storing information about the last timestep and acceleration 
for each particle, we have not implemented any of this. 



2.4. Time Step 

The time-step is subjected to the following constraints: In accord with the 
Courant-Friedrichs-Lewy (CFL) condition for stability of finite difference meth- 
ods [35], we require 

At = Chydio M & l(\ut+c s y (30) 

where Chydro < 1 is the CFL number, and m and c s are the fluid velocity in the ith 
direction and sound speed of the flow respectively. In presence of a source term, S, 
we modify the estimate of the fluid velocity according to 

\Sj\Ax 

W + Cs ~~ * [(N + c s Y + 2 \Si\ Ax]V2 _ + Cs ) 

where Si is the component of the source term affecting the velocity u,. For the 
purpose of accuracy, rather than stability, for the collisionless particles we likewise 
require 

„ a(t) Ax 
Af = C part W 32 
Max(|«j|) 

with Cp ar t < 1 and with £/, corrected as in Eq. (31) but with c s = 0. Finally, 
we require that the background expansion remains limited during each integration 
cycle. This allows us to time center the value of a(t) in our integration schemes 
above and neglect its changes with time. We thus enforce 



At < C cxp - 
a 



(33) 



with C cxp ~ (1 - 2) x 10" 
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2.5. Code Units 

The natural choice for the dimensional units of the above physical equations is 
given by the following lengths, mass and time scales 

L* = L box (34) 

p* = p c ^l m (35) 

U = Hq 1 (36) 

whereLbox is the size of the computational box and p c = 3Hq/8ttG = 1.879 x 
10~ 29 h 2 is the critical density of the universe, with h = 7?o/100 km s _1 Mpc _1 . 
The units for the other quantities are defined in terms of these as 

u* = Hq 1 L box (37) 

P* = p*u\ (38) 

4>* = < (39) 

T* = m proton P^/ksP* (40) 

where m pro ton is the proton mass and fee Boltzmann's constant. 



3. ADAPTIVE MESH REFINEMENT APPROACH 

In this section we outline the structure of the AMR algorithm. After introducing 
the formal notation, we describe in some detail the scheme for hydrodynamics with 
self-gravity and its modifications when a collisionlcss component is also included, 
in the simple case of two levels of refinement. We then describe the extension of 
the scheme to the general multilevel case. 

3.1. Multilevel Volume Discretization, Variables and Operators 

The underlying discretization of the D-dimensional space is given as points 
(io, ...,Zd-i) = i £ The problem domain is discretized using a grid T c lP 
that is a bounded subset of the lattice. T is used to represent a finite-volume dis- 
cretization of the continuous spatial domain into a collection of control volumes: 
i 6 T represents a region of space, 

Vi = [ih,{i + u)h] (41) 

where h is the mesh spacing, and u 6 TlP is the vector whose components are all 
equal to one. We can also define face-centered discretizations of space based on 
those control volumes: T e = {i ± ^e d : i E T}, where e d is the unit vector in the 
d direction. T ed is the discrete set that indexes the faces of the cells in T whose 
normals are e d : 

A i+led = [(t + e d )h, (i + u)h],i + \e d e T ed . (42) 
2 Z 

We define cell-centered discrete variables on T: 



4> : r -> K ; 
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and denote by <f>i E W n the value of <t> at cell i e T. We can also define face-centered 
vector fields on T: 

F=(F ,...,F O - 1 ) ,F d :T ed ^R m 
and define a discretized divergence operator on such a vector field: 

1 D_1 

{D ■ F)i = - ^2 (F d , iHe * F d>i _, ed ),i e T. (43) 

We will find it useful to define a number of operators on points and subsets of IP . 
We define a coarsening operator by: C r : IP — > Z D , 



Cr(t) = 



( 






id-i 


) 




r 


, . . . , 


r 





where r is a positive integer. These operators acting on subsets of Z can be 
extended in a natural way to the face-centered sets: C r (T e ) = (C r (T)) e . For any 
set T C T, we define G(T, r), r > 0, to be the set of all points within a | • |-distance 
r of T that are still contained in T: 

G(T,r) = rn U T + i 

\i\<r 

where \i\ — max (\id\)- We can extend the definition to the case r < : 

g(T,r)=T-G(T-T,-r). 

Thus Q(T,r) consists of all of the points in T that are within a distance — r from 
points in the complement of T in T. In the case that there are periodic boundary 
conditions in one or more of the coordinate directions, we think of the various sets 
appearing here and in what follows as consisting of the set combined with all of 
its periodic images for the purpose of defining set operations and computing ghost 
cell values. For example, Q(T,r) is obtained by growing the union of T with its 
periodic images, and performing the intersections and differences with the union of 
r with its periodic images. 

We use a finite-volume discretization of space to represent a nested hierarchy 
of grids that discretize the same continuous spatial domain. We assume that our 
problem domain can be discretized by a nested hierarchy of grids T° ...T lmax , with 
r' +1 = C^, 1 (T l ), and that the mesh spacings h l associated with T l satisfy jj+y = 

n l re j. The integer n l re j is the refinement ratio between level I and I + 1. These 
conditions imply that the underlying continuous spatial domains defined by the 
control volumes are all identical. In this paper we will further assume n l re j is even. 
In the case where there are only two levels, we will refer to them as coarse and fine, 
with the notation {I = 0, 1 = 1} — > {c, /}, and n° e j — ► n re f. 

We make two assumptions about the nesting of grids at successive levels. We 
require the control volume corresponding to a cell in O i 1 is either completely 
contained in the control volumes defined by il l or its intersection has zero volume. 
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We also assume that there is at least one layer of level-/ cells separating level-(Z + 1) 
cells from level-(7 — 1) cells: Q(C n i ^(f^ +1 ), 1) C ft 1 . We refer to grid hierarchies 
that meet these two conditions as being properly nested. 

From a formal numerical analysis standpoint, a solution on an adaptive mesh hi- 
erarchy {^'};™5 X approximates the exact solution to the Partial Differential Equa- 
tions only on those cells that are not covered by a grid at a finer level. We define 
the valid region of SI 1 as 

n l valid = n l -c nl jn l+1 ). 

A composite array ip comp is a collection of discrete values defined on the valid 
regions at each of the levels of refinement: 

We can also define valid regions and composite arrays for face-centered variables: 
^vaiid = — C n i ^(Q l+1 ' ed ). Thus, Qyaiid consists of d-faces that are not cov- 
ered by the d-faces at the next finer level. A composite vector field F comp = 
^pl.validylmax is defined as follows: 

Thus a composite vector field has values at level I on all of the faces not covered 
by faces at the next finer level. 

We want to define a composite divergence D comp {F l+1 > comp ,F l ' comp )i for i G 
^vaiid- To do this, we construct an extension of F l ' comp to the edges adjacent 
to £l l va i id that are covered by fine level faces. On the valid coarse-level cZ-faces, 
^d = F d comp . On the faces adjacent to cells in Sl l vatid , but not in fi^ id , we set 
F l d = (F l d +1 ' comp ), the average of F l d +1 onto the next coarser level: 

\*d )ii + ke d - ( n ,)D-1 2^ b d,i+\e* > ll + 2 E (d,+ U Q- ■ 
v reJ ' i+ie«ef 

Here T d is the set of all fine level d-faces that are covered by A i[+ i e d. Q± consists 
of all the rf-faces in fl l on the boundary of with valid cells on the low (± = — ) 
or high (± = +) side: 

C^ 1 - {i ± \e d : i ± e d e ti valid , i e C^JQ^ 1 )}. 
Given that extension, our composite divergence is defined as: 

D comp(pl+l,comp^ pl,compj. = D . p! , i £ il l vaHd . (44) 

It is useful to express D comp as the application of the level divergence operator 
D applied to extensions of F l ' comp to the entire level, followed by a step that 
corrects the cells in $l l valid that are adjacent to Q l+1 . We define a flux register 
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5F l+1 associated with the fine level 

5F l+1 = (SFl ) +\...,SF}+\) 

Let F l be any coarse level vector field that extends F l > com v i. e . 



pi _ Tjl,comp o'< ed 



Then, for i E n l vaHd , 

D comp(j?l+l,comp t pLcomp^, = + D R (SF l+1 )i, (45) 

where SF l+1 is a flux register set to be 

^m = (^+ 1 )-^onC^ + UCl_ 

and Dr is the reflux divergence operator, given by the following for valid coarse 
level cells adjacent to ft' +1 : 

(2=0 ±=+,-: 

For the remaining cells in ftj, aiid , Dr(5F 1+1 ) is defined to be identically zero. 

We can use this notation to define the discretizations of Poisson's equation we will 
be using to compute self-gravity. On a single level, ft', we define A', the discrete 
Laplacian, to be the standard 2D + 1 point operator, with the values used on ghost 
cells computed using quadratic interpolation: 

Ay = D ■ F 1 ' 4, (46) 

Fd ~ — u — (47) 

4 = T{$\$- X )i for i e dn l (48) 

where <9ft' = £/(ft', 1) — ft • Here the interpolation function T is an 0(h 3 ) estimate 
of the value on the ghost cell obtained from interpolating from values of (j) 1 ^ 1 on 
^valid ano - from the values of 4> l on ft' ; for details, see [36] . We can then define the 
composite Laplacian A comp applied to all of the valid data on all levels, in terms 
of that operator and refluxing operations. 

(A™'^ = (Ay^ 4 ); + D R (SF 1 ^) for i e Q l valid (49) 
§ pl,4> = (pl+ht} _ pl,4> (50) 

where (fy l ' ext is some extension of (f> l < com P to all of ft'. The resulting operator depends 
only on the valid values of </> in the grid hierarchy (modulo roundoff considerations; 
cf. Ref. [37] for further details). 
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3.2. AMR for Compressible Flows with Self- Gravity 

The starting point for this work is the algorithm described in [2] for solving 
hyperbolic conservation laws on nested refined grids. For the case of two levels, we 
assume that the solution on both levels is known at time t c . The basic steps to 
evolve the solution on both levels to time t c + At c can be summarized as follows: 

1. Update the solution on the coarse grid: 

U c (t c + At c ) = U c (t c ) - At c (D ■ F c ) + At c S(U c ) on n c . (51) 

Here S(U C ) is computed as in Eq. (15) (with the body force, /, set to zero), and 
the discrete fluxes F are local functions of U c . We also initialize flux registers 
associated with using the same fluxes 

SF f = -F c . 

2. Advance the solution from to +At f on the fine grid n re f times, n re fAt^ = 
At c : 

U f (t f + At f ) = U f (t f ) - At f {D ■ F f ) + AtS(U f ) on Q f (52) 
SF f += —(Ff) 

7l re f 

t f += At f . 

Any values required to compute the stencil that are contained in Tf — Qf are 
computed by interpolating the coarse grid values U c (t c ), U c (t c + At c ), using linear 
interpolation in time, and piecewise linear interpolation in space. 

3. Synchronize the values at the old and new times: 

U c (t + At c ) = (U f (t + At c )) on C nref (n f ) (53) 
U c (t + At c )+= D R (5F) 

where {■)% denotes the arithmetic average onto the coarse cell i of all of the values 
defined on fine grid cells contained in i. 

To extend this algorithm to the case of self-gravity, we must solve the Poisson's 
equation for the gravitational potential due to the mass distribution of the fluid 
on the coarse and fine levels. As usual the coarse level is advanced first, and the 
solution at t c + At c is used to provide time interpolated boundary conditions for 
the fine level at intermediate time steps. In the case of hyperbolic equations the 
finite characteristic speeds ensure that a fully consistent multilevel solution can 
be recovered at synchronization time with the refluxing operation. Because of its 
elliptic nature, however, in order to preserve its multilevel character Poisson's equa- 
tion should be solved simultaneously on all levels. Notice that the coupling among 
levels is enforced by the continuity of the potential (Dirichlct) and of its normal 
derivative (Neumann) at the coarse/fine grid interface. Therefore, to maintain the 
multilevel character of Poisson's equation when the levels are not synchronized yet, 
we obtain a single level solution to Poisson's equation on the coarse level and apply 
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a lagged estimate of the effect of the coarse/fine matching conditions at refinement 
boundaries, following the ideas developed in [25, 26, 38] for incompressible fluids. 
This leads to the following modifications to the algorithm given above. 

0. At simulation start, when all levels are synchronized, we compute a composite 
grid solution to Poisson's equation 

(L corop -V comp )(t c )i = P c (t c )i i e n c valid (54) 

{L comp ' f ^ comp )(t c )i = p f (t c )i ieCl f (55) 
as well as the coarse grid solution 

{L c <t> c ){t c ) = p c (t c ) on n c . (56) 

Here, and in what follows, we will denote by L(t) = j^A with superscripts I, c, 
/, comp indicating the particular discretization of the Laplacian operator. 

1. Together with U(t c ), the acceleration f c (t) derived from the composite so- 
lution of the potential on the coarse grid is used to compute the coarse grid fluxes 
and source terms. After updating the conserved quantities using (51), we compute 
the coarse grid potential at the new time 

(L c <j) c )(t c + At c ) = p c (t c + At c ) (57) 

where in the latter step we have approximately corrected the coarse grid single-level- 
solution of the potential for the effects due to the solution on the finer grid [25, 26]. 
We use the solution in Eq. (58) with (j) c ' com ' p (t c ) to obtain boundary conditions 
interpolated in time for the potential at the fine level at intermediate timesteps. 

2. We apply the update (51) on the finer level. At tf = t c in order to compute 
the hydrodynamic fluxes and source terms we use the force, f*, derived from the 
composite potential solution to Eq. (55). At intermediate steps, t c < tf < t c + 
Ai c , we solve the following Poisson's equation on the fine grid with interpolated 
boundary conditions 

(L f <j) f )(t f ) = p f (t f ) on Q f (59) 
(l> f (t f )=I[<p f (t f ),4> c {t f )] on 8Sl f 

where <j) c {t^) is obtained by linear interpolation in time as 

4> c (t f ) = {l-a) 4> c > comp {t c ) +a4> c > comp (t c + At c ) , a = 

We compute the forces at the new timestep and apply the correction (17)-(18) to 
the fluid momentum and kinetic energy. 

3. At time of synchronization, = t c + At c , we solve the set of equations 
(54)-(56), for the composite and single level solution of the potential. We use the 
composite solution to derive the new force, f n+ , and apply momentum and kinetic 
energy corrections (17)-(18) to obtain time centered forces on both coarse and fine 
levels. The flow of the calculation restarts from step 1 with the gravitational force 
known at all levels. 



! Please write \titlerunninghead{< (Shortened) Article Title>} in file ! 15 



3.3. AMR with Particles 

Due to the time refinement character of the AMR technique the solution on 
different levels is advanced with different timesteps. This implies that the density 
field represented by the particles evolved on the finer level may not be available 
on the coarser level unless they are synchronized. However, this is information is 
necessary to solve Poisson's equation. Therefore, we find it convenient to introduce 
effective particles to recover such information in a computationally inexpensive way 
and without compromising the code accuracy. 

Thus we introduce an aggregation operation V — > (V) 1 that projects a collection 
of particles covered by fl l onto a set of effective particles, with no more that one 
particle per cell. If p E (V) 1 , then 

m p = ^ m p i (60) 
p'.x p , eVi 

x p = — nip'Xpi (61) 

mp P ':x p ,eVi 

u p = Yl m p ,u p' ( 62 ) 

The aggregation operation conserves the monopole and dipole terms but causes 
information to be lost on the quadrupole moment of the aggregated particles [39] , 
which provides corrections to the potential of order h 2 . Thus the aggregation 
step preserves second order accuracy. Note, also, that the potential and force 
fields obtained through the aggregated particles are only used to provide boundary 
conditions for the finer level. 

Restricting again the discussion to the case of two levels of refinement, the changes 
in the algorithm described above are given as follows: 



0. At simulation start, we partition the particles into ones that will be evolved 
using the coarse and fine time steps. If V is the set of all particles, 

V s = {p G V : Xp G G(Sl f , -n ref n buf )} (63) 
V c = V -V s . (64) 



The parameter n&„/ is chosen so that the support for the Particle-Mesh interpo- 
lation function, used to calculate the force acting on the particle and the particle 
mass distribution on the grid, is completely contained in ft? for all of the fine grid 
time steps, n re fAt^ — At c . Using n^uf = 1 and C par t — 0.5 is sufficient for all 
choices of Particle-Mesh scheme used here (see Appendix 1). 



We then define the set (7 3 ^) c of fine particles aggregated on the coarse grid using 
Eq. (60)-(62). 
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Finally, in computing the gravitational potential the densities in Eq. (54)- (56) are 
modified to account for the mass distribution of the particles: 



P w) = P r md in + E ™A {i + l z u)h : c - x v {tc) 



p{{t c ) = p{' nui \t c )+ E m p w\ 
+ E m p wi 




hf 

(i+\u) hf -x p {t c ) 
hP 



(65) 



(66) 



Here W is one of the Particle-Mesh assignment schemes used to spread the particle 
mass on the grid, described in Appendix 1. Note that the addition to p c of the 
density field due to {V^Y has no effect on the composite solution, because the 
support of W for each particle <E (P^) is contained entirely in C„ rc/ (iV). 

1. The gravitational force resulting from the composite solution of the potential 
is used to compute both the fluid fluxes and source terms, as well as to perform 
the update (22)-(23) for the particles in V c U {V^Y- For the latter, we interpolate 
the accelerations from the grid cells to particle positions with one of the methods 
described in Appendix 1. 

2. At the end of each fine time step, while + At* < t c + At c , we compute the 
new fine potential, (fif , modifying the mass density as 



P f( t f) i = pf^ id (tf) i + e ™ P w( {i+1 * u)h ^ f Xp{tf) ) 

P evf V / 

V- ljr f{i+ \u)h f -x p {t f )\ . . 

+ E m P w ( - — ~\c — — , « e 



(67) 



where the positions of the particles in V c at the intermediate times are given by 
linear interpolation between x p (t c ) and x p (t c + At c ). We then use the acceleration 
due to this field to update the fine particle velocities using (24), and the fine fluid 
state using (17)-(18). 

3. The synchronization step is analogous to the one in Sec. 3.2: we calculate a 
single grid and composite grid solution of the potential using the total mass density 
distribution of fluid and particles given by Eq. (65)-(66). The gravitational force 
derived from the composite potential is used to apply the corrections to the particle 
velocity, given in Eq. (24), and to the fluid momentum and kinetic energy, given 
in Eq. (17)-(18), at all levels. Finally, the sets V c , V f and (T f Y are upgraded 
following the definitions in (63)- (64) to account for the new particle positions. The 
flow of the calculation restarts from step 1. 



3.4. The General Multilevel Algorithm 

For the case of a general (Z mox + l)-level calculation, we assume that at simulation 
start all the particles have been partitioned into groups corresponding to the levels 
on which they shall be advanced, and that the particles being advanced by finer 
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grids have been aggregated into effective particles for the next coarser levels. This 
is summarized in the procedure: 
for / = 0, . . . l m ax - 1 do 

V l+1 = {xeP l : x p {t l ) e g(Q t+1 ,-n ref n buf )} 

p l <_ -pi _ -pi+i 

end for 

■pUa* = 

for / = l max - 1, ... ,0 do 

v l c = (r l+1 up l c +i y 

end for 

We can then describe the algorithm for advance (I) that advances the solution 
at level I, < I < l max by one time step. We assume that the solution is known at 
time t l + At 1 in a process that includes a recursive application of advance. At the 
beginning of advance, we assume that: we know the fluid state U l and particle 
state P l associated with that level at time t l ; we know the fluid, particle, and 
potential at the next coarse level at times t 1 ' 1 , t 1 ' 1 + At'" 1 ; a composite solution 
as well as single level solution of the gravitational potential at time t l has been 
computed on level I — 1 and on all finer levels; 

We then advance the solution by a time step At 1 , with n^jAt 1 — At 1 ' 1 , as 
follows: 

1. Using U(t l ) along with the accelerations, f\, i £ Q l , computed from 4> l , we 
calculate the hyperbolic fluxes F l and the source terms S(t l ). We update the 
conserved quantities 

U(t l + At 1 ) = U(t l ) - At l D ■ F l + At l S l (t l ) (68) 
6F i+i = _pi if l < lmax 

6F l += (F l ) if I > 

and advance the positions and velocities of the particles in V 1 U V l c using (22) and 
(23). 

2. With t l <— t l + At 1 , we solve Poisson's equation on level I 

(L l 4> l )(t l )i = P nmd '\t\ + £ m p w( {i + ^ U)h ^ Xv[tl) \ (69) 



p£p l 



^ / (i + \u)h l - xJt l )\ ^ / (i + \u)h l -xJt l )\ , 

+ J2m p w[ { 2 ; £Li + ^ m P W( { - 2 > ,ien l 



p€p l c \ / pGp' 



cl> l (t l ) i = i[<t> l (t l ),4> l - 1 (t%,i£dn l 

where </> (f) is obtained by linear interpolation in time between 0' _1 > com P(t' _1 ) 
and ft-hcompfjl-i + A^" 1 ) with the latter defined in Eq. (58). As in Eq. (58) we 
then estimate 

4> l > com P(t l ) = <f> l {t l ) + (^■ com P(i'- 1 ) - cp^t 1 ' 1 )) (70) 
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which is used to provide boundary conditions for the potential at the next finer 
level, 4> l+1 {t). 

3. We call advance recursively n l re j times. 
t l+1 = t l 

while t l+1 < t l + At 1 do 

advance (I + 1) 

t i+i = t i+i + At l+i 

end while 

4. If t'- 1 < t l < t^+At 1 ' 1 , we obtain {<t> l '' comp (t l )}i>>i by solving the composite 
equations : 

{L c ° m *' 1 ' 4> comp W)i = P fluid ' l '(t l )i + ]T m p w( {i+llU)H ^ Xp{tl) \ (71) 

pEV 1 ' V / 

_ / (i + \u)h 1 ' -x p (t l ) \ 

+ L m p w [ — —f^l 

* ^ ^valid i 1=1, l m ax 

<p l ' comp (t l ) i =l[(j> l ' comp (t l ),4> l - 1 (t%, i e dn l . 

The field so obtained is used to compute accelerations at the new time, and update 
the fluid and particle velocities using (17) - (18) and (24) respectively. 

5. The solution at level I is synchronized with the solutions at the finer levels. 

U(t l + At 1 ) = (U(t l+1 )) onC n , m i+1 ) 

re J 

U{t l + At l )+= D R (5F l+1 ) 

We upgrade the sets V 1 ' for V = l — l, ■■■,l m ax an d (V l ' +1 ) 1 ' for I' = I — 1, l max ~ 1, 
according to the new particle positions. The flow of the calculation restarts from 
step 1. 

4. CONVERGENCE TESTS 

We have implemented the above schemes in a Cosmological Hydromagnetic AMR 
Radiation Many-body (CHARM) code. The code is based on the CHOMBO AMR 
library and it is implemented in a hybrid CH — [-/Fortran77 language. Additional 
physics modules, such as radiation [40], cosmic-rays [41, 42] and magnetohydrody- 
namics will be presented elsewhere. In the following, we focus on numerical tests 
to assess the performance of the algorithms in terms of accuracy and applicability 
to problems of direct interest. Performance tests will be presented elsewhere. 

Unless explicitly stated otherwise, in the following we use these CFL coefficients 
for the time step: Chydro = C part = 0.5 and C oxp = 0.01. In addition, we restrict 
the results to the case of a TSC interpolation scheme which, in accord with previous 
authors, we find to give the most accurate results. 
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TABLE 1 
Convergence tests: collisionless case 



-/Vpart 


Li 


Ri 


L 2 


i?2 


Loo 


Roo 








position 








8 


1.3e-07 


1.9 


1.4e-07 


1.9 


1.9e-07 


1.8 


16 


3.5e-08 


2.0 


3.9e-08 


2.0 


5.4e-08 


2.0 


32 


8.9e-09 


2.0 


9.8e-09 


2.0 


1.3e-08 


2.1 


64 


2.2e-09 


2.0 


2.4e-09 


1.9 


3.0e-09 


1.6 


128 


5.6e-10 




6.2e-10 




9.8e-10 










velocity 








8 


1.3e-04 


1.9 


1.4e-04 


1.9 


1.9e-04 


1.8 


16 


3.5e-05 


2.0 


3.8e-05 


2.0 


5.3e-05 


2.0 


32 


8.8e-06 


2.0 


9.7e-06 


2.0 


1.3e-05 


2.1 


64 


2.2e-06 


2.0 


2.4e-06 


1.9 


3.0e-06 


1.6 


128 


5.5e-07 




6.2e-07 




9.7e-07 










force 








8 


6.5e-02 


1.9 


7.1e-02 


1.9 


9.4e-02 


1.8 


16 


1.7e-02 


2.0 


1.9e-02 


2.0 


2.7e-02 


2.0 


32 


4.4e-03 


2.0 


4.8e-03 


2.0 


6.5o-03 


2.1 


64 


l.le-03 


2.0 


1.2e-03 


1.9 


1.5e-03 


1.6 


128 


2.8e-04 




3.1e-04 




4.8e-04 





t We use equal number of cells and particles, iVpart = AT ce n, a cell-centered force scheme, and 
a constant A£ = 1.6 X 10~ 4 . 



Errors and convergence rates are calculated as follows. At a given resolution, r, 
for any given cell or particle, i, we estimate the error on a computed quantity, 
with respect to the analytic solution,^ (i), as 

5q r (i) = q c r (i)-q a r (i). (72) 

We then compute the n-norm of the error, i.e. 

l/n 



L n (5q r ) = \\Sq r \\n = []T \Sq r (i)\ n Vi\ " (73) 



where v\ is either the i-th cell volume or the inverse of the number of particles; 
finally we estimate the convergence rate as 

_ ln[L n (Sq r )/L n (Sq s )} 

"~ MAXr/AXs) ■ 1 j 

For the cases studied below we report the L l7 L 2 and norm of the errors. 



4.1. Zel'dovich's Pancake 

We begin with a classical test problem for cosmological codes, the evolution 
of a one-dimensional plane-wave perturbation in an expanding background. In 
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FIG. 1. L2 norm of the error in position (top), velocity (center) and force (bottom) as 

a function of the number of grid cells. Left panels correspond to the case in which the number 
of particles and grid cells is the same, N ce u = N var t- Central and right panels correspond to 
N P art = \N ce u, with particles initially placed either at cell nodes (middle panels) or at cell 
centers (right panels). In all cases a two point cell centered force stencil is used. See legend for 
the meaning of the symbols (UG= Uniform Grid, a is the expansion parameter). 



Zel'dovich's formulation [43], the comoving position and peculiar velocities of col- 
lisionless matter evolve as 

X{t) = q+ % p{q) (75) 
v(t) = a(t) x(t) (76) 



where q,p are the Lagrangian initial position and displacements, a(t) is the ex- 
pansion factor and b(t)/a(t) describes the growth factor of the perturbation. For a 
closed universe (0 TO = 1), a(t) = (3H a t/2) 2 ? 3 (Eq. [2]) and b(t) = 2a 2 /5 [43]. 
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TABLE 2 

Convergence tests: collisionless, variable timestep, linear phase 



A^part 


Li 


Ri 


L 2 


i?2 


Loo 


Roo 








position 








8 


6.4e-06 


1.9 


6.9e-06 


1.9 


9.2e-06 


1.8 


16 


1.7e-06 


2.0 


1.9e-06 


2.0 


2.6e-06 


2.0 


32 


4.3e-07 


2.0 


4.7c-07 


2.0 


6.4e-07 


2.1 


64 


l.le-07 


2.0 


l.le-07 


1.9 


1.4e-07 


1.6 


128 


2.5e-08 




2.8e-08 




4.5e-08 










velocity 








8 


8.9e-04 


1.9 


9.6e-04 


1.9 


1.3e-03 


1.8 


16 


2.3e-04 


2.0 


2.6e-04 


2.0 


3.6e-04 


2.0 


32 


5.9e-05 


2.0 


6.5e-05 


2.0 


8.9e-05 


2.1 


64 


1.5e-06 


2.0 


1.6e-05 


1.9 


2.0e-05 


1.6 


128 


3.7e-07 




4.2e-06 




6.6e-06 










force 








8 


6.5e-02 


1.9 


7.1e-02 


1.9 


9.4e-02 


1.8 


16 


1.7e-02 


2.0 


1.9e-02 


2.0 


2.7e-02 


2.0 


32 


4.4e-03 


2.0 


4.8e-03 


2.0 


6.5o-03 


2.1 


64 


l.le-03 


2.0 


1.2e-03 


1.9 


1.5o-03 


1.6 


128 


2.8e-04 




3.1e-04 




5.1e-04 





t We use equal number of cells and particles, -/V pa rt = JVcelli a cell-centered force scheme, a 
variable timestep, ^| = C eX p(%), C e x P = 10~ 2 , and a = 0.0221. 

Setting the initial displacements to a sinusoidal form, p(q) = 5Asm(kq)/2, where 
k is the perturbation wavenumber, we obtain 

x(t) = q + aAsm(kq) (77) 
v(t) = a a Asin(kq) (78) 
p(t) = po [1 + a A k cos(kq)}^ 1 . (79) 

The solution described by Eq. (77) becomes singular, that is dx/dq = when 
^collapse = (Ah) -1 . At this stage, particles trajectories cross at x = q = n/k 
and a caustic forms. In the following we use a sta rt = 1/51, a C oiiapse = 1/2 and 
k = 2ir/h~ 1 L hox . 

4- 1.1. Collisionless Component 

Table 1 demonstrates the second order accuracy of the code. The different 
columns report, as a function of numerical resolution, the L\, L 2 and norms 
of the error on the particles position, velocity and force, and the corresponding 
convergence rates of the errors. The L 2 norm of these errors is also reported graph- 
ically in the left hand side panels of Fig. 1. For this case we use equal number 
of cells and particles, N part — N ce n, a cell-centered force scheme, and a constant 
= 1.6 x 1CP 4 . In Fig. 1 we also report the L 2 norm of the errors for the case 
in which = ^N cc n and the particles are initially placed either at cell nodes 
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TABLE 3 

Convergence tests: collisionless, variable timestep, nonlinear phase 







uniform 


grid 








AMR 






iVpart 


Li 


-Ri 




Roo 


Li 


Ri 




— R 


: 

I'Tnax 










position 










8 


3.0e-02 


1.7 


4.6e-02 


1.6 


2.9e-02 


1.9 


4.5e-02 


1.9 


l 


16 


9.5e-03 


1.8 


1.6e-02 


1.8 


7.9e-03 


2.0 


1.2e-02 


2.0 


l 


32 


2.7e-03 


1.9 


4.7e-03 


1.8 


2.0e-03 


2.0 


3.1e-03 


1.9 


2 


64 


7.2e-04 


2.0 


1.3e-03 


2.0 


5.1e-04 


2.0 


8.3e-03 


1.9 


2 


128 


1.8e-04 


2.0 


3.3e-04 


1.8 


1.2e-04 


2.0 


2.2e-04 


1.3 


3 


256 


4.6e-05 




9.2e-05 




3.0e-05 




9.2e-05 




3 










velocity 










8 


5.9e-02 


1.4 


9.4e-02 


1.2 


5.1e-02 


1.8 


7.9e-02 


1.9 


1 


16 


2.2e-02 


1.6 


4.1e-02 


1.2 


1.4e-02 


1.8 


2.1e-02 


1.5 


1 


32 


7.4e-03 


1.6 


1.8e-02 


1.3 


4.1e-03 


2.1 


7.4e-03 


2.2 


2 


64 


2.4e-03 


1.7 


7.3e-03 


1.4 


9.6e-03 


2.1 


1.6e-03 


1.9 


2 


128 


7.4e-04 


1.8 


2.8e-03 


1.6 


2.2e-04 


2.0 


4.5e-04 


1.3 


3 


256 


2.1e-04 




9.1e-04 




5.7e-05 




1.8e-04 




3 



force 



8 


1.5e-01 


1.1 


2.6e-01 


0.7 


l.le-01 


1.3 


1.8e-01 


0.9 


1 


16 


7.4e-02 


1.0 


1.6e-01 


0.5 


4.5e-02 


2.4 


9.3e-02 


2.6 


1 


32 


3.7e-02 


1.0 


1.2e-01 


0.6 


8.2e-02 


1.2 


1.6e-02 


0.1 


2 


64 


1.8e-02 


1.2 


7.8e-02 


0.7 


3.5e-03 


0.5 


1.5e-02 


2.3 


2 


128 


8.1e-03 


1.4 


4.6e-02 


1.0 


2.5e-03 


1.6 


3.0e-02 


1.3 


3 


256 


3.1e-03 




2.2e-02 




8.2e-04 




1.2e-02 




3 



t We use equal number of cells and particles, -/V pa rt = JVcelli a cell-centered force scheme, a 
variable timestep, ^| = C exp (j), Cexp = 1CT 2 , and a = 0.479. 



(central panels, NCP for node centered particle) or cell centers (right panels, CCP 
for cell centered particle). Both these configurations lead to a slower convergence 
rate than our reference case illustrated in the left hand side panels. 

In cosmological simulations, the timestep during the initial stages is determined 
by the expansion rate of the background. Thus in Table 2 we report the same 
quantities as in Table 1 but for the case of a (variable) timestep determined by 
s| = Cexp(f ), with C exp = 10~ 2 . The L 2 norm of these errors are also reported 
as filled symbols in the left hand side panels of Fig. 2. The errors were computed 
after ten timesteps so that the system is still in the linear regime. And, in fact, the 
convergence rates are the same as in Table 1. The same correspondence in terms 
of convergence rates also exists for the case in which N part = |iV ce ii, as illustrated 
for the NCP case by the star symbols in right hand side panels of Fig. 2. 

Next, we consider the errors during the nonlinear regime of the calculation. In 
particular, we consider the solution just prior to the caustic formation, when the 
background expanded by a factor 25 since the simulation start. On the left hand 
side of Table 3 we report the L\ and error norms and convergence rates as in 
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Table 2 while the L — 2 errors are shown as open symbols in the left hand side panels 
of Fig. 2. We see that in the nonlinear regime the convergence rates of the errors in 
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FIG. 2. L2 norm of the error in position (top), velocity (center) and force (bottom) as a 
function of the number of grid cells. Left panels correspond to the case in which the number of 



particles and grid cells is the same, N ce n 



part ■ 



Right panels correspond to N par 
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with particles initially placed at cell nodes. See legend for the meaning of the symbols (UG= 
Uniform Grid, a is the expansion parameter). 



the particle positions, velocities and forces have worsened in a minor, appreciable 
and considerable way, respectively. 

Finally, we test the performance of the AMR code. We use a constant refinement 
ratio, n ro f = 2, and refine cells enclosing a mass larger than 1.5 the average value. 
A maximum of three levels of refinement were allowed. All runs used a first level of 
refinement for about 30% of the calculation, except for the lowest resolution case 
for which the percentage was 12%. The second level of refinement was only used by 
the three higher resolution runs and only for about 5% of the time. Finally, a third 
level of refinement was employed only by the two finest runs and only for less than 
1% of the time. Similarly, finer grids cover a progressively smaller fraction of the 
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computational domain. The L\ and L m error norms and convergence rates for the 
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FIG. 3. Comparison between two numerical experiments, employing a uniform grid (left) 
and two levels of refinement (right), in terms of error in the particle position (top), velocity 
(middle) and force (bottom). The initial set up includes 32 grid cells (bound by vertical lines) and 
32 particles uniformly distributed. Only one quarter of the grid is shown, focusing on the critical 
region where the caustic forms. Vertical lines indicate cells' boundaries. 



AMR runs are reported on the right hand side of Table 3 while the L — 2 errors are 
shown as spur symbols in the left hand side panels of Fig. 2. These results show 
that employing AMR during the nonlinear evolution improves the convergence rate 
of the solution in such a way that they resemble the values in the linear stage. This 
is true in this example for the errors in the position and velocity and to a lesser 
extent, the force, which is more affected by coarse/fine boundary effects. 

Fig. 3 compares the errors in position, velocity and force for a fixed grid (left) 
and an AMR grid (right) calculations. It focuses on the region where, and the 
times when, the caustic forms and AMR operates. Thus only one quarter of the 
total grid is shown (with the cell boundaries of the base grid indicated by vertical 
lines and marked by integer labels), and the various errors are plotted only for the 
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last one third of the simulation run. At the beginning of that time span (when 
a ~ 0.2), the particles have clustered sufficiently at the grid center and a first level 
of refinement is created. Six grid meshes are refined (only three between boundaries, 
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FIG. 4. Same as Fig. 1 but for the case in which a staggered scheme is used for the force 
calculation. 



12-16, are shown in Fig. 3). The second level of refinement is created much later 
and only affects one base (or two refined) grid mesh(es) for the last ten per cent of 
the simulation. Fig. 3 shows that the errors in position, velocity and force of the 
particles close to the point where the caustic forms are much reduced when AMR 
is employed. The generation of a level of refinement is accompanied by a change in 
the mass distribution and the potential field. When this happens, a particle may 
experience a sudden change in terms of the force field acting upon it. These effects 
are responsible for the somewhat 'errant' behavior of the force error, as illustrated 
in the bottom-left panel of Fig. 3. Overall, though, for each particle the force 
fluctuations in the AMR case are significantly smaller than the force errors in the 
uniform grid case. 
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As a last result for this section, in Fig. 4 we plot L 2 convergence errors analogous 
to Fig. 1 but for the case in which the force was computed with a staggered scheme. 
Comparison of the two figures shows that when N palt — N cc u the solution obtained 
with the staggered scheme also converges with second order accuracy, while being 
characterized by smaller errors. This is in agreement with previous findings [44]. 
However, when the number of particles is halved (-/Vp art = ^N ce n) the results ob- 
tained with a staggered force scheme worsen more dramatically than for the cell 
centered case, showing a very poor convergence rate. Closer inspection shows that 
when the particles are sparse on the grid, oscillations appear in the potential due to 
the discrete character of the matter distribution as reproduced on the grid. Since 
this affects only the quality of the staggered scheme the problem may be related 
to the inconsistency of this scheme with the centering of the stencil used for the 
discretization of the Laplacian operator. We shall return to this issue in the next 
test case where the problem reappears with more dramatic effects. 

4-1.2. Collisional Component 

We now turn to the performance of the hydrodynamic part of the code. We use 
the same tests employed in the previous section for the the collisionless component. 
However, since in this case the solution refers to an Eulerian grid, in order to have 
the velocity and density at a given grid location x(q, t) from Eq. (78)- (79) we must 
invert Eq. (77). We first consider the errors in the linear regime, in analogy 



TABLE 4 

Convergence tests: Godunov's scheme* 



iVpart 


Li 


Ri 


L 2 


R 2 


Loo 


Roo 








density 








8 


2.7e-05 


2.2 


2.9e-05 


1.9 


4.2e-05 


1.2 


16 


5.7e-06 


2.2 


7.6e-06 


2.2 


1.8e-05 


1.8 


32 


1.2e-06 


2.0 


1.7e-06 


2.2 


5.3e-06 


1.9 


64 


3.0e-07 


2.0 


3.8e-07 


2.1 


1.4e-06 


1.9 


128 


7.4e-08 




8.8c-08 




3.7e-07 










velocity 








8 


3.3e-05 


2.0 


3.6e-05 


2.0 


5.2e-05 


2.0 


16 


8.1e-06 


2.0 


9.0e-06 


2.0 


1.3e-05 


2.0 


32 


2.0e-06 


2.0 


2.2e-06 


2.0 


3.2e-06 


2.0 


64 


5.0e-07 


2.0 


5.6e-07 


2.0 


7.9e-07 


2.0 


128 


1.2c-07 




1.4e-07 




2.0e-07 










force 








8 


1.6e-02 


2.0 


1.7e-02 


2.0 


2.4e-02 


2.0 


16 


3.9e-03 


2.0 


4.3e-03 


2.0 


6.1e-03 


2.0 


32 


9.7e-04 


2.0 


l.le-03 


2.0 


1.5c- 03 


2.0 


64 


2.4e-04 


2.0 


2.7e-04 


2.0 


3.8e-04 


2.0 


128 


6.0e-05 




6.7e-05 




9.5e-05 





t We use a cell-centered force scheme and a constant timestep,-^ = 1.6 X 10 4 . 
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TABLE 5 

Convergence tests: Godunov's scheme, linear regime 



iVpart Ll Rl L2 R2 Loo Roc 



density 

8 2.0e-04 2.2 2.1e-04 1.9 3.0e-04 1.1 

16 4.2e-05 2.1 5.7e-05 2.1 1.4e-04 1.7 

32 9.5e-06 2.0 1.3e-05 1.9 4.3e-05 1.4 

64 2.4e-06 1.9 3.5e-06 1.5 1.6e-05 0.9 

128 6.5c-07 1.3e-06 8.4e-06 



velocity 

8 2.1e-04 2.0 2.3e-04 2.0 3.3e-04 2.0 

16 5.1e-05 2.0 5.7e-05 2.0 8.2e-05 2.0 

32 1.3e-05 2.0 1.4e-05 2.0 2.0e-05 2.1 

64 3.2e-06 2.0 3.5e-06 2.0 5.0e-06 1.9 

128 7.9e-07 8.9c-07 1.3e-06 



force 

8 1.5e-02 2.0 1.7e-02 2.1 2.3e-02 2.0 

16 3.6e-03 2.0 4.0e-03 2.0 5.7e-03 2.0 

32 9.0e-04 2.0 1.0e-04 2.0 1.4e-03 2.0 

64 2.2e-04 2.0 2.5e-04 2.0 3.5e-04 2.0 

128 5.6e-05 6.2e-05 8.9e-05 



t We use a cell-centered force scheme, a variable timestep, ^ = C e xp(j )> C eX p = 10 2 , and 
a = 0.0221. 



to Table 1 and Table 2 for the collisionlcss component. In Table 4 we report the 
Li, L 2 and L^ norms of the error for the gas density, velocity and force and the 
corresponding convergence rates, for the case of fixed time step, |j = 1.6 x 10~ 4 . 
Similarly, the left hand side of Table 5 reports the L\, L 2 and L x errors for the 
case in which the time step is set by the background expansion rate. The L 2 errors 
for the fixed and varying timestep cases arc also shown by the filled symbols in the 
left and right hand side panels of Fig. 5, respectively. These tests show the second 
order accuracy of the implemented scheme. We note, however, that in the case 
of variable time step the convergence rate of the L^ norm of the density error is 
slower. Next, the left hand side of Table 6 reports the errors and convergence 

rates for the case of a uniform grid calculation, well in the nonlinear regime, close 
the formation of the caustic (a = 0.479). As before, the L 2 errors are also shown as 
open symbols in Fig. 5. Unlike the collisionlcss case, here is the convergence rate 
of the density that is mostly affected, particularly at low resolutions (cf. [31, 45]). 
As illustrated by the L^, norm, the error is dominated by the contribution of a few 
cells, located where the caustic forms. 

Finally, we test the performance of the AMR code. As for the collisionless com- 
ponent, we use a constant refinement ratio, n re f = 2, refine cells enclosing a mass 
larger than 1.5 the average value and allow for a max of three levels of refinement. 
The results at a = 0.476 are reported in the right hand side columns of Table 6 for 
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TABLE 6 

Convergence tests: Godunov's scheme, nonlinear regime 









unifnrm 










AMR 






iVpart 


Li 


Ri 


Loo 


Roo 


Li 


Ri 


Loo 


Roo 


Imax 










density 










8 


2.2e-01 


0.5 


8.2e-01 




1.6e-01 


0.7 


8.2e-01 


— 


l 


16 


1.6e-01 


0.7 


1.2e00 




1.0e-01 


0.6 


1.5e00 


— 


l 


32 


9.7e-02 


0.9 


1.5e00 


0.0 


6.4e-02 


2.0 


2.4e00 


1.3 


2 


64 


5 Op-02 


1.6 


1.5e00 


1.0 


1.6e-02 


2.0 


l.OeOO 


1.1 


2 


128 


1.6e-02 




7.3e-01 




3.9e-03 




4.7e-01 




3 










velocity 










8 


2.7e-02 


1.4 


9.4e-02 


0.4 


1.6e-02 


1.8 


8.6e-02 


0.7 


1 


16 


1.0e-02 


1.4 


7.0e-02 


0.4 


4.7e-03 


1.7 


5.3e-02 


0.6 


1 


32 


3.8e-03 


1.5 


5.3e-02 


0.6 


1.4e-03 


2.4 


3.6e-02 


2.2 


2 


64 


1.3e-03 


1.6 


3.6e-02 


0.9 


2.6e-04 


1.8 


9.6e-03 


1.9 


2 


128 


4.3e-04 




2.0e-02 




7.2e-05 




3.1e-03 




3 










force 










8 


6.8e-02 


1.6 


2.3e-01 


0.5 


3.7e-02 


1.7 


2.0e-01 


0.9 


1 


16 


2.3e-02 


1.5 


1.6e-01 


0.5 


l.le-02 


2.3 


l.le-01 


0.8 


1 


32 


8.4e-03 


1.5 


l.le-01 


0.6 


2.2e-03 


2.0 


6.3e-02 


1.7 


2 


64 


2.9e-03 


1.6 


7.3e-02 


0.9 


5.5e-04 


1.8 


1.9e-02 


2.2 


2 


128 


9.5e-04 




3.9e-02 




1.6e-04 




4.0e-03 




3 



t We use a cell-centered force scheme, a variable timestep, ^ = C eX p(^), Cexp = 10 2 , and 
a = 0.479. 

the Li, errors and convergence rates. There we also show the maximum level 
employed by each run. L 2 errors are reported as spur symbols in the right panels 
of Fig. 5. 

The use of refined grids in terms of fraction of grid covered and fraction of the 
simulation time is very similar to the corresponding collisionless case. Similar is 
also the benefit of AMR, which improves the convergence of the solution to rates 
very similar to those characterizing the linear regime. This is indeed a powerful 
performance of the AMR technique. 

4.2. Effect of C exp on the Solution Quality 

We have investigated how the error depends on the choice of the parameter C exp 
for the above problem. In particular we have computed the error accumulated 
during an interval Aa <C a co u apse , for values of the expansion parameter a = 0.196 
and a — 0.091, and for values of C exp ranging from 10~ 3 to 0.5. We consider both 
the fluid and the collisionless case. We use a uniform grid with 32 zones on a side 
and, for the collisionless case, we use one particle per cell. The results are rather 
independent of the norm type. We find that the errors introduced in the particles 
velocity and position is quite stable, except for the largest values of C exp . On the 
other hand, the errors in the fluid components decrease steadily as C exp is reduced, 
spanning a factor ~ 5 before reaching a plateau for C exp < (1 — 2) x 10~ 2 . 
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FIG. 5. L2 norm of the error in density (top), velocity (center) and force (bottom) as a 
function of the number of grid cells. See legend for the meaning of the symbols (UG= Uniform 
Grid, a is the expansion parameter). 



4.3. Homologous Dust Cloud Spherical Collapse 

In this section we test the ability of the code to follow the collapse of a pressure- 
less (dust) sphere of matter [46]. The problem is described by the following equa- 
tions 

^ =-G^± (80) 



dt 2 



M 



(81) 

M 



with initial conditions 



P^- Q > - \ if r(i)>R (82) 

u(i,i = 0) = (83) 
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2 4 2 4 

r/ Ax r/ Ax 

FIG. 6. Phase space distribution of particles for the spherical collapse of a pressureless 

cloud. Cyan open circles correspond to the analytic solution and filled triangles to the numerical 
simulation result. The right panel is for the case of a uniform grid whereas the left panel correspond 
to an AMR calculation with two levels of refinement. 



where M(r) is the mass enclosed within a distance r from the sphere center, R is 
the radius of the sphere and f(r) is a function (with f'(r) < 0) that depends solely 
on r. For the hydrodynamic case null density will be approximated with a value 
p(r > R) <C min[p(r)]. The problem admits a self similar solution in implicit form 
which reads 

(1 - £) 1/2 e /2 + sin-^l - 1/2 = t (84) 
r = £r , u=— (£ -1) , r = t\ J (85) 

where (p) r — 3Af(r)/47rr 3 is the average density within a radius r. 

From the numerical point of view, the problem is challenging in two respects: 
during collapse the force potential becomes progressively steeper and, therefore, 
more demanding for the gravity solver. In addition, since the problem has inherent 
radial symmetry and we are solving it on a Cartesian grid, the ability of the code 
at preserving that symmetry will be tested. 

For the collisionlcss component we initially set particles with null velocity at 
the center of cells whose distance from the cloud center is less than R (=1). This 
produces a homogeneous density distribution everywhere inside R, except close to 
the cloud edge due to the discreteness of the grid. In Fig. 6 we compare the position 
of each particle in phase-space (v r ,r) as given by the code (black filled triangles) 
with the analytic solution (cyan open circle). The left panel corresponds to the 
case in which a uniform grid is used whereas for the right panel solution AMR was 
employed. We allow for two levels of refinement and tag cells with a total mass four 
times as high as the initial value. The comparison between analytic and numerical 
solution in Fig. 6 is made for a number of evolution times, expressed in terms of the 
adimensional collapse time t co u = it/2. Noticeably, the code follows the particles 
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FIG. 7. Phase space distribution of particles for the spherical collapse of a pressureless 

cloud. The initial cell-to-particle ratio is 4. Cyan open circles correspond to the analytic solution 
and filled triangles to the numerical simulation result using a cell-centered (left) and face-centered 
(right) force scheme, respectively. 



motions with high accuracy all the way down to the time of collapse. In particular, 
there is no sign of artificial asymmetries. Additional levels of refinement were 
dynamically generated towards the final phase of the collapse. At the latest time 
shown (r = 0.98 t co /;), we can see the improvement due to the employment of finer 
grids in the collapsing region. Note that towards the edge of the cloud the particles 
are trailing. This is due to the inability to reproduce a perfectly homogeneous 
sphere near the cloud edge from the beginning. The region affected by this is about 
one mesh size wide. 

When exploring the accuracy of cell-centered versus face-centered force schemes 
our tests suggest that, again, in the uniform grid case the latter perform slightly 
better, at the level of ca 15%. In analogy with the analysis of Sec. 4.1, we have 
tested this further, for the case in which the ratio of cells to particles is significantly 
larger than one. This situation may easily occur depending on the adopted criterion 
for refinement and on the efficiency for the generation of the refined grid out of the 
tagged cells. In Fig. 7 we compare the solutions obtained with a cell-centered (left) 
and a staggered (right) force scheme for an initial cell-to-particle ratio of 4. The 
initial dust sphere is placed on a grid of 256 cells on a side, and 8 particles are 
aligned along its radius out to 32 cells from its center. The code output is plotted 
at three different solution times close to the time of collapse t — t co u. This test 
shows that when the number of cell-to-particle ratio is significantly higher than 
one, the staggered force scheme tends to produce spurious results. This is in line 
with the findings in the previous test problem in section 4.1. On the other hand, 
the cell-centered scheme seems well behaved. We find that the qualitative result 
does not change when we use a two or four point force stencil, when the number of 
particles to resolve the sphere is changed or when the sphere center is shifted by a 
fraction of a mesh size in an arbitrary direction. 
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Thus the staggered force scheme, although apparently more accurate than its 
cell-centered counterpart when the number of cells is comparable to, or less than, 
the number of particles [44], it gives rise to spurious results when particles are 
sparse on the grid. Therefore, caution must be exercised when employing staggered 
schemes for force evaluation. 
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FIG. 8. Density and velocity profile for spherical collapse of a pressureless cloud of 

gas. Cyan open circles correspond to the analytic solution and black filled dots to the numerical 
simulation results. 

Finally, the results for the collisional case are illustrated in Fig. 8. The plot 
compares the density (top) and velocity (bottom) profiles of the numerical solution 
(threaded black dots) and the analytic solution (cyan dots). The latter extend only 
out the cloud size, whereas the numerical solution includes the region covered by the 
finest level. Two times during the collapse are shown: t — 0.72t co u and t = 0.92r co /;, 
corresponding to the low and high curves, respectively. At these times one and two 
levels of refinement have been generated, respectively. The chosen times are close 
to the collapse time, when the errors have accumulated and the simulation becomes 
more challenging. Nevertheless, as for the collisionless component, the code follows 
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accurately the evolution of the density and velocity profiles of the collapsing cloud. 
Again, close to the cloud edge the density profile is smoother and the velocity field 
slower than the analytic solution. The size of the region affected by this is again of 
order of the coarse mesh size and it is partially ascribed to the crude representation 
of the cloud edge on the grid. 

4.4. Energy Conservation: Layzer-Irvine Equation 

For pure hydrodynamics conservation of the total (kinetic+thermal) energy is 
enforced by our conservative Godunov's method. When gravity is added, energy 
conservation should still hold, but is not explicitly enforced in our scheme. Finally, 
with an expanding background energy is not conserved. For a collection of parti- 
cles that interact only gravitationally, the evolution of the energy of the system is 
regulated by the Layzer-Irvine equation, which reads 

±[a(t)(£ + W)] = -a£ (86) 

where £ is the kinetic energy associated to the peculiar motions of the particles 
and W their gravitational potential energy due to the overdensity produced by 
their mass distribution. Clearly in absence of expansion (a = 1, a = 0) Eq. (86) 
reduces to the ordinary energy conservation equation. Otherwise it describes the 
change in the total energy of the system due to the adiabatic expansion of the back- 
ground. The derivation and physical meaning of Eq. (86) is reviewed in Ref. [27]. Its 
applicability to hydrodynamic simulations is discussed in, e.g., Ref. [32], in which 
case a monoatomic gas is assumed and £ includes both the kinetic and thermal 
energy of the gas. 

Eq. (86) can be integrated in time giving 

a[£(a) + W(a)} - a [£(a ) + W(a )] = - f £da. (87) 

J ao 

We can evaluate the integral on the RHS of the above equation numerically with 
the trapezoidal rule and assess the accuracy of the code at tracking the energy of 
the system through the quantity 

_ a[£(a) + W(a)} - a [£(a ) + W(a )} + /° £da 
£ [a W(a ) - aW(a)] ' [ ' 

We first test the energy conservation accuracy of the code for the case of the 
collapse of a pressureless cloud. This is the problem studied in the previous section. 
We carry out three AMR calculation with different base grid sizes, namely 16, 32, 
64 corresponding to 4, 8, 16 cells per cloud radius respectively. In these runs, cells 
enclosing more then four times the initial mass content were tagged for refinement 
and a maximum of two 2 refinement levels were allowed. The results of the test 
are reported in the left hand side panel of Fig. 9 where we plot the error in the 
total energy, Se, as a function of resolution, for both the particle (open) and the 
gasdynamic (filled) case. The plots show that with 16 zones per cloud radius the 
error in the energy is at the level of a per cent or so. 
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FIG. 9. Energy conservation error based on Eq.(88) for the collapse of a pressureless 

spherical cloud (left) and a cosmological run (right). See text for details. 



Next we test code accuracy at tracking the energy of the system in a cosmological 
run. For the purpose we use a A-Cold Dark Matter cosmology with parameters 
fl m =0-3, Oa = 0.7, fib = 0.04, for the energy density in total matter, dark energy 
and baryonic matter respectively; and H = 70 km s _1 Mpc -1 for the Hubble 
constant. The physical domain has a size of L=91.43 Mpc on a side. We execute 
three runs with different numerical resolution. The first two runs employ a uniform 
grid with 32 3 and 64 3 cells, respectively, and the same number of particles as grid 
cells. The third run uses a base grid with 32 3 cells and 32 3 particles, and two 
additional levels of refinement created in region where the total mass enclosed in 
a cell exceeds the initial value by a factor cigth. The initial conditions where 
generated on a 64 3 grid and coarse-averaged to a 32 3 grid for the low-resolution- 
uniform and AMR runs. 

The results are presented in the right panel of Fig. 9 where we plot Se (top), £ 
(middle) and W (bottom) as a function of expansion parameter a for each resolution 
case. The plots show that when using a uniform grid the total energy of the system 
is evolved with an accuracy at the percent level (~ 2% and 1% for the 32 3 (dot) 
and 64 3 (dash) cases, respectively), with most of the error generated at startup. In 
the AMR case (solid line), however, our error parameter Se increases visibly when 
refinement levels are created (at a ~ 0.15 and a ~ 0.2 for the first and second level 
respectively). The reason for this is simple. When a level of refinement is created 
the potential energy of the system changes suddenly (see solid and dot lines in the 
bottom panel) throwing off the balance bewtween particle/gas velocities and their 
potential energy. As a result a large error in the sense of Eq.(88) is generated. This 
is so, even though with the additional level of refinement the potential energy of the 
system is more accurate as it gets closer the value from the high resolution run (see 
solid and dash lines in the bottom panel). Over time the kinetic energy readjusts to 
the new potential (middle panel; notice that the internal energy is negligible) and 
a balance between the two forms of energy is reestablished. Because the potential 
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energy associated with the newly formed structures is larger than that of the system 
at the time when the refinement levels were first generated, the new balance between 
kinetic and potential energy reduces substantially the error as time progresses. At 
simulation end the AMR run (with a base grid of 32 3 cells) produces estimates of 
5e, £ and W very close to the high resolution run. 

4.5. Santa Barbara Galaxy Cluster 

In this section we carry out the calculation defined by the 'Santa Barbara Cluster 
Comparison Project' [47] and compare the results of our code with those from 
different codes implemented independently by other authors and based either on 
similar or different techniques. 

The problem consists of simulating the formation of a galaxy cluster in a Standard 
Cold Dark Matter universe. The cosmological parameters assumed were fl m = 1 an d 
fib = 0.1 for the total and baryonic mean mass density in units of the critical density, 
respectively; H n — 50 km s _1 Mpc -1 for the Hubble constant; as — 0.9 for the 
present-day linear rms mass fluctuation in spherical top hat spheres of radius 16 
Mpc; for the baryon density. The computational domain has a size of L=64 Mpc on 
a side. The initial matter fluctuation are characterized by a power spectrum with 
an asymptotic spectral index, n = 1, and are 'constrained' so that at simulation 
end a massive structure has formed at the center of the computational box. 

The simulation was initialized at z = 40 with two grids already in place: a base 
grid covering the entire 64 Mpc 3 domain with 64 3 cells and 64 3 particles; and a 
second grid, also with 64 3 cells and 64 3 particles, but only 32 Mpc on a side and 
placed in the central region of the base grid, thus yielding an initial cell size of 
0.5 Mpc. Refinement is applied only in this central, higher resolution region and 
is based on a local density criterion: cells with a total mass of 6.4 x 10 10 M Q or 
more were refined. We allowed for 5 levels of refinement (for a total of a 6 levels 
hierarchy), with a constant refinement ratio n re f — 2. The size of the finest mesh 
is about 15 comoving kpc. We use the following CFL coefficients for the time step: 
Chydro = 0.8, C part = 0.5 and C oxp = 0.02. 

At simulation end a halo finder based on the spherical overdensity method [48] 
was run in order to define the center of the galaxy cluster. The radial profiles for six 
quantities of interest are presented in Fig. 10 together with results from two other 
simulation codes: ENZO, which is an Eulerian AMR code similar to ours, and HYDRA 
(as run by Jenkins & Pearce), which combines smoothed particle hydrodynamics 
(SPH) and adaptive particle-particle-particle-mesh (AP 3 M) method for the N-body 
part [49]. These two codes are meant to be representative of the high resolution 
grid based and SPH approaches, respectively. 

Results are shown down to scales of about 30 kpc which is just above the nominal 
resolution at finest level of refinement. The plot shows that there is good agree- 
ment among the results of the different codes particularly with ENZO (even though 
we used one less level of refinement and use a different refinement criterion). The 
discrepancies, when significant, are consistent with those already found in the ex- 
tensive comparison paper in Ref. [47]. In fact, there is very good agreement in 
terms of dark matter density distribution (top left) which is well fit by the analyt- 
ical form proposed in Ref. [50] with parameters specified in the caption. Similarly, 
there is good matching of the solution in terms of gas density distribution, except 
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FIG. 10. Radial profile of dark matter (top left), baryonic gas (top right) temperature 
(middle left), baryonic fraction (middle right), radial velocity for dark matter (bottom left) and 
gas (bottom right). In addition to the results from CHARM (open dots), for comparison we also 
show those from the ENZ0 AMR code (filled triangles) [21] as well as those from the HYDRA SPH 
code (open stars) [49]. 



in the inner regions within 100 kpc, where the two AMR solutions flattens and the 
SPH solution keeps on increasing. More significant is the difference in temperature 
distribution (middle left), which drops in the SPH solution for the inner regions 
and stays constant for the AMR case. These differences were already found and 
discussed in Ref. [47]. (See also Ref. [19] for similar findings.) Their origin is not 
fully clear, although as suggested in Ref. [47], it could be ascribed to the different 
way in which SPH and grid based methods treat shocks. 

Next panel (middle right) shows the profile for the ratio of gas to dark matter 
mass, normalized to the global value. While our solution is in good agreement with 
ENZO's and deviates from HYDRA's in the inner regions and is somewhat in between 
the two beyond 1 Mpc or so. There is significant scatter in the results from the full 
set of codes found in Ref. [47], at the level of 0.1-0.2. Nevertheless, it is pointed 
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out in Rcf. [47] that when integrated out to the virial radius of the system (2.7 
Mpc), the SPH codes seem to predict a systematically slightly smaller values for 
this quantity than high resolution grid based codes. The reason for this is still not 
clear. 

Finally, both the gas and dark matter radial velocity profiles agree quite well 
at all radii. Some differences may arise due to slight differences in the simulation 
timing, as pointed out in Ref. [47]. Note that at larger radii (last few points), 
typically characterized by wider scatter, both ENZO's and HYDRA's results tend 
to be below the average value defined in Rcf. [47] . 

5. CONCLUSIONS 

We have presented a new code based on AMR technique for systems comprising 
collisional and collisionless components coupled through a long range force. We 
have thus extended the scheme in [2] to include collisionless particle dynamics 
and gravity arising from the mass distribution of the two components. For the 
hydrodynamics we use a slightly modified dircctionally unsplit Godunov's method 
based on Ref. [28] . As for the collisionless component we have implemented various 
time centered modified symplectic schemes based on both the kick-drift-kick and 
drift-kick-drift sequence. Our implementations of these schemes appear to perform 
comparably. We have also used several types of stencil to calculate the force from 
the potential. We find that while the staggered schemes appear more accurate 
when the number (density) of particles is at least as large as the number of grid 
cells, it produces spurious results when the particles are sparse on the grid. Cell 
centered stencils thus seem more reliable, especially when a five point cell centered 
discretization of the Laplacian operator is used. 

Due to the time refinement character of the AMR technique the solution on dif- 
ferent levels is advanced with different timesteps. Synchronization issues then arise 
as the multilevel solution to the elliptic equation needs to be solved simultaneously 
on all levels. In particular, the density field represented by the particles evolved 
on finer levels may not be available on coarser levels when they are not synchro- 
nized. Similarly, one cannot account for the effects of the mass distribution on 
finer levels on the multilevel solution of the potential, unless all levels are synchro- 
nized. Among other features of the code, we have thus introduced an aggregation 
procedure to cost-effectively represent on the coarser levels the particles on finer 
levels without compromising the code accuracy and performance. We have also 
implemented a procedure for estimating (when the coarse and finer levels are not 
synchronized) the effects on the coarse potential produced by the matching con- 
ditions at refinement boundaries between coarse and fine solution, as they would 
arise in a full multilevel calculation. We performed several standard tests which 
illustrate the code accuracy as well as the advantages of the AMR technique for the 
study of both self gravitating hyperbolic systems, collisionless system and hybrid 
systems. 
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APPENDIX: CHARGE ASSIGNMENT SCHEMES 

Indicating with Ax the mesh size in one dimension we find the nearest grid point 
scheme (NGP) defined for order r = 1 as 

ttt ( ~ x f 1 if \x-x\ < Ax 12 .. 

W(x - x) = i , 1 . 1 ~ ' (A.l) 

v I otherwise v ; 

in which the assigned charge distribution is discontinuous as the particle cross the 
cell boundary; the cloud in cell (CIC) scheme defined for r = 2 as 

nrt- \ [l-JHZil if \x-x\<Ax . . n . 

W(x - x) = < Ax 1 1 ~ (A.2) 

I otherwise 

in which the assigned charge distribution is continuous but the first derivative is 
not; the cell boundary; and the triangular shape cloud (TSC) scheme defined for 
r = 3 as 



W(x - x) 



, , 2 

3 / \x— x 



4 



if \x-x\< Ax/2 
i_(|_M_^ 2 if Ax/2 < \x - x\ < 3Ax/2 ( A ' 3 ) 

otherwise 



in which both the assigned charge distribution and first derivative are continuous. 
These schemes retain their properties when they are extended to a multidimensional 
case in the form of a product 

D 

W(Z-x) = l[W j (x j -x j ). (A.4) 

3=1 
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